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Cholesky  factorization  that  is  also  called  the  Shur  or  Bcrlekamp-Massey 
algorithm  . 

We  review  how  fast  Cholesky  factorizations  are  usually  used  to  compute  re 
flection  coefficients  from  correlations.  Then  we  show  how  the  factorizations 
may  be  run  backwards  to  compute  Cholesky  factors  and  correlations  from 
reflection  coefficients.  This  generalizes  a  result  usually  attributed  to 
Robinson  and  Treitel.  One  of  our  main  purposes  is  to  emphasize  that  the 
Levinson  recursions  for  going  back  and  forth  between  correlations,  reflection 
coefficients,  and  autoregressive  filter  parameters  may  be  replaced  with  a 
dual  set  of  recursions  for  going  back  and  forth  between  correlations, 
reflection  coefficients,  and  moving  average  filter  parameters. 

When  the  correlation  matrix  R  describes  an . autoregressive  moving  average 
time  series,  then  each  column  of  the  Cholesky  factor  of  R  consists  of  a 
p-dimensional  Kalman  gain  vector  in  its  first  p  non-zero  entries,  followed 
by  a  homogeneous  recursions  for  subsequent  values.  This  means  the  fast 
Cholesky  factorization  of  R  may  be  further  speeded  up  by  using  these  Kalman 
gains  to  set  initial  conditions  in  a  homogeneous  recursion.  It  means  also 
that  the  Kalman  gains  obey  the  same  vector  recursions  as  the  columns  of  the 
Cholesky  factor  of  R.  This  fast  Kalman  gain  algorithm  may  be  read  out  of  the 
LeRoux-Gueguen  algorithm  or  out  of  the  Friedlander  algorithm. 
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the  Cholesky  factor  of  R.  This  fast  Kalaan  gain  algoritha  nay  be  read  out 
of  the  LeRonx-Gnegnen  algoritha  or  out  of  the  Friedlander  algoritha. 

1 .  Introdnction 

The  Levinson  recursions  provide  an  efficient  algoritha  for  factoring 
the  inverse  of  a  Toeplitz  correlation  aatrix  R  into  its  upper  and  lover 
triangular  Cholesky  factors.  These  factors  produce  a  Graa-Schaidt 
orthogonal ization  of  the  underlying  tiae  series.  The  recursions  are 
routinely  used  to  coapute  reflection  coefficients  for  iapleaenting  whitening 
and  predicting  filters  in  lattice  fora.  Thvy  are  also  used  to  go  back  and 
forth  between  reflection  coefficients,  order-increasing  whiteners,  and 
correlations.  In  this  paper  we  tell  the  dual  of  this  story,  based  on  a 
factorization  of  the  correlation  aatrix  R,  itself.  One  of  our  aain  purposes 
is  to  show  how  the  Levinson  recursion  for  going  back  and  forth  between 
correlations,  reflection  coefficients,  and  order- increasing  whiteners  aay  be 
replaced  with  a  dual  set  of  recursions  for  going  back  and  forth  between 
correlations,  reflection  coefficients,  and  order  increasing  synthesizers.  ' 

Ve  use  the  Levinson  recursions  to  derive  the  LeRoux-Gueguen  algoritha 
[1]  for  factoring  a  Toeplitz  correlation  aatrix  R,  colunn-by-coluan.  We 
organize  the  algoritha  into  a  coupled  set  of  vector  recursions  and  show  that 
it  is,  indeed,  a  lattice  algoritha.  By  re-ordering  the  colnan  variables 
into  rows,  we  convert  the  LeRoux-Gueguen  algoritha  into  Friedlander 's 
lattice  algoritha  [2]  for  factoring  a  Toeplitz  correlation  aatrix,  row-by¬ 
row.  (Conversely,  we  could  have  re-ordered  Friedlander’s  algoritha  to 
obtain  the  LeRoux-Gueguen  algoritha.)  This  discussion  shows  the  two 
algorithas  to  be  two  different  ways  of  looking  at  a  fast  Cholesky 
factorization  that  is  also  called  the  Shur  [3],  Berlekaap-Nassey  [4],  or 


Abstract 


Te  use  the  Levinson  recursions  to  derive  the  LeRoux-Gueguen  algorithm 
for  factoring  a  Toeplitz  correlation  matrix  R,  coluam-by- column .  Ve 
organize  the  algorithm  into  a  coupled  set  of  vector  recursions  and  show 
that  it  is,  indeed,  a  lattice  algorithm.  By  re-ordering  the  column 
variables  into  rows,  we  convert  the  LeRoux-Gueguen  algorithm  into 
Friedlander ' s  lattice  algorithm  for  factoring  a  Toeplitz  correlation  matrix, 
row-by-row.  This  shows  the  two  algorithms  to  be  different  ways  of  looking 
at  a  fast  Cholesky  factorization  that  is  also  called  the  Shur  or  Berlekaap- 
Massey  algorithm. 

Ve  review  how  fast  Cholesky  factorizations  are  usually  used  to  compute 
reflection  coefficients  from  correlations.  Then  we  show  how  the 
factorizations  may  be  run  backwards  to  compute  Cholesky  factors  and 
correlations  from  reflection  coefficients.  This  generalizes  a  result 
usually  attributed  to  Robinson  and  Treitel.  One  of  our  main  purposes  is  to 
emphasize  that  the  Levinson  recursions  for  going  back  and  forth  between 
correlations,  reflection  coefficients,  and  autoregressive  filter  parameters 
may  be  replaced  with  a  dual  set  of  recursions  for  going  back  and  forth 
between  correlations,  reflection  coefficients,  and  moving  average  filter 
parameters. 

When  the  correlation  matrix  R  describes  an  autoregressive  awving 
average  time  series,  then  each  column  of  the  Cholesky  factor  of  R  consists 
of  a  p-dimensional  Kalman  gain  vector  in  its  first  p  non-zero  entries, 
followed  by  a  homogeneous  recursion  for  subsequent  values.  This  means  the 
fast  Cholesky  factorization  of  R  may  be  further  speeded  up  by  using  these 
Caiman  gains  to  set  initial  conditions  in  a  homogeneous  recursion.  It  means 
also  that  the  Kalman  gains  obey  the  same  vector  recursions  as  the  columns  of 


Morf  [5]  algorithm.  When  we  rnn  the  recursion  backwards  to  obtain  Cholesky 
factors  and  correlations  from  reflection  coefficients  we  generalize  the 
result  of  Robinson  and  Treitel  [6]  for  obtaining  correlations  from 
reflection  coefficients. 

When  the  correlation  matrix  R  describes  an  autoregressive  moving 

average  (ARHA)  time  series,  then  each  column  of  the  Cholesky  factor  for  R 

consists  of  a  p-dimensional  Kalman  gain  vector  in  its  first  p  non-zero 

entries,  followed  by  a  homogeneous  recursion  for  subsequent  values  [7]. 
This  means  the  fast  Cholesky  factorization  of  R  may  be  further  speeded  up  by 
using  these  Kalman  gains,  together  with  ARKA  parameters,  to  set  initial 
conditions  and  run  a  homogeneous  recursion  to  generate  columns  [7].  It 

means  also  that  the  Kalman  gains  obey  the  same  vector  recursions  as  the 
columns,  themselves  [7],  [8].  Thus  the  fast  Kalman  gains  may  be  read  out  of 
the  LeRoux-Gueguen  algorithm,  or  out  of  the  Friedlander  algorithm,  by 

reading  internal  variables  out  of  a  lattice.  The  algorithm  may  be  slightly 
modified  to  produce  the  Morf,  Sidhu,  Kailath  algorithm  [9]  for  computing 
Kalman  gains  on  a  fixed  length,  time-varying  lattice  [11]. 

The  diagram  in  Figure  1  summarizes  our  perspective  and  our  results. 

The  right  side  of  the  diagram  we  take  to  be  well  understood:  Levinson 

recursions  may  be  used  to  factor  R-*  and  to  go  back  and  forth  between 

correlations  {rn},  order  increasing  whiteners  {a11},  prediction  error 

2 

variances  {on  }  and  reflection  coefficients  {knJ .  The  left  size  of  the 
diagram  we  take  to  be  less  well  understood  and  it  is  our  intention  to 
clarify  it  and  extend  known  results  concerning  it.  To  this  end  we  show  how 
the  LeRoux-Gueguen  and  Friedlander  algorithms  may  be  used  to  factor  R  and  to 
go  back  and  forth  between  correlations  (rn),  order  increasing  synthesizers 
{ha} ,  prediction  error  variances  {oa^},  and  reflection  coefficients  {ka). 


Faster  Cholesky 
factorization 


Fast  Kalman  Gains 


Le  Roux-Gueguen 

< - 

R  =  HD2H' 


Levinson 

- > 

-1  -2 
R  =  A'D  A 


Figure  1. 


Connections  between  Algorithms 


These  results  are  used  to  speed  up  the  Cholesky  factorization  in  the  ARMA 
case  and  to  read  the  fast  Kalman  gain  algorithm  out  of  the  factorization. 


2.  Linear  Statistical  Models  for  Stationary  Sequences 
Let  x  =  (x0,  xj, . . . '  denote  t  real  random  variables  drawn  from  a 
wide  sense  stationary  random  sequence  with  mean  value  sequence  zero  and 
correlation  sequence  (rt).  The  first  two  moments  of  x  are 


Ex  =  0 


Exx'  =  R  «  {r . .  . 


|i-i| 


r0  rl 


t-1 


r0  rl*  * ‘ rt-2 


ri  ro 


Call  R  a  (txt)  Toeplitz  correlation  matrix. 
Define  the  exchange  matrix  J: 


J  = 


0 

1 


1 


0 


=  J* 


JJ  =  I 

This  matrix  turns  column  vectors  up-side  down  and  row  vectors  right-side 
left : 

Jx  =  . . . i  x0) 

X  J  -=  (  Xj.j  #  ....  Xe ) 

The  vector  x  is  drawn  from  a  wide-sense  stationary  sequence  that  does  not 


know  forvard  time  from  reverse  time.  Therefore  Jx  has  the  same  moments  as 


This  says  R  is  sot  only  Toeplitz,  but  also  centro-  or  J-symmetric. 


2.1  Synthesis 

The  correlated  vector  x  may  be  synthesized  from  an  nncorrelated  vector 

Da: 

x  =  HDu  ;  Eon*  =  I 
R  =  HD2H' 

The  diagonal  matrix  D  and  triangular  matrix  H  are  defined  as  follows: 


Ve  call  HD2H*  a  lower-upper  (LU)  Cholesky  factorization  of  R  and  x  =  HDu  a 
forward  synthesis  of  x. 

The  super  and  subscripting  of  elements  in  H  is  somewhat  arbitrary,  but 
generally,  hj  tells  how  the  random  variable  u^  is  used  to  build  the  random 
variable  Xj+j* 

The  vector  Jx  may  be  synthesized  as  follows: 


Jx  -  JHJJDJJu 


E  =  GC2G* 


The  matrices  G  and  C  are  defined  as  follows: 


C2  =  JDJ  =  diag  (o2_lt  ....  a2,  o2) 


G  =  JHJ  = 


Ve  call  GC^G'  an  npper-lower  (UL)  Cholesky  factorization  of  R  and  Jx  = 
a  backward  synthesis  of  Jx. 


2.2  Analysis 

The  random  vector  x  may  also  be  analyzed  to  produce  the  nncorrelated 
random  variables  Du  in  a  Gram-Schmidt  procedure: 

Ax  *  Du 
ARA’  =  D2 

The  triangular  matrix  A'  is  defined  as  follows: 


1 


1 

a. 


2 

a. 


a  backward  analysis  of  Jx.  The  first  of  these  interpretations  comes  from 
the  following  manipulation  of  BRB'  =  C^: 


(BRB')"1  =  C~2 
R-1  =  B'C"2B 


2.3  Somnary 

From  these  factorizations  of  R  and  R-*  we  have  these  farther 
connections  between  analysis  and  synthesis  parameters: 


H  =  A  1 

A'  =  H  1  * 

G  =  B-1 

B'  =  G"1' 

The  LU  and  UL  factorizations  of  R  may  therefore  be  rearranged  as  follows: 

RA'  =  HD2 

RB'  =  GC2 

These  resnlts  for  analysis  and  synthesis,  and  for  factorization  of  R 
and  R-1,  are  summarized  in  Fignre  2.  It  is  already  obvious,  and  onr  fast 
algorithms  will  confirm  it,  that  H(G)  and  B'(A')  play  dual  roles.  This 
leads  us  to  suspect  that  recursions  for  columns  of  B'(A')  should  lead  to 


similar  recursions  for  columns  of  H(G) 


Synthesis 


Analysis 


Forward 

x  =  HDu  Ax  =  Du 

R  =  HD2H'  (LU)  R_1  *  A'D-2A  (DL) 

RA'  =  HD2  R_1H  =  A'D-2 

Backward 

Jx  =  GJn  BJx  =  CJu 

R  =  GC^1  (DL)  R-1  -  B'C-2B  (LD) 

RB'  =  GC2  R_1G  =  B'C-2 

Connections 

H  =  JGJ  (L)  A'  «  JB'J  (U) 

G  *  JHJ  (D)  B'  -  JA'J  (L) 

D  =  JCJ  D_1  =  JC_1J 

C  =  JDJ  C”1  *  JD-1J 

H  =  A'1  A'  =  H-1' 

G  =  B"1  B'  =  G"1' 

Figure  2.  Forward  and  Backward  Analysis  and  Synthesis 


Begin  with  the  correlation  matrix  R  and  the  following  rearrangements  of 
its  UL  and  LU  factors: 

RA'  =  HD2 
RB'  =  GC2 

3.1  Levinson  Recnrsions 

Call  Jan+2  the  vector  of  non-zero  elements  in  the  (n+1)  column  of  A* 
and  an+*  the  vector  of  non-zero  elements  in  the  (n+1)  column  of  B',  counting 
from  the  right-most  column  of  B'.  Then,  from  the  factorization  above,  we 
have  the  relations 


Note  what  happens  when  Jan+1  and  an+1  are  replaced  by  the  previous 


columns  Ja11  and  a11: 


n+1 


L  i  n+l-i 
i=0 


o  1  n+1 

2 

a 

n 

• 

rl 

• 

n 

a 

0 

• 

• 

• 

• 

0 

a 

r  >4  •  •  •  1 4  r 

n+1  1  o  J 

0 

»  « 

n 

Cn+1  . 

Choose  Jan+1  end  an+1  to  be  the  following  coupled  recursions  in  order  to 
satisfy  the  factorizations: 


0 

1 

_  n+1 
Ua  J 

as 

J*n. 

+  kn+l 

***•1  o. 

n+1 

n 

0  ' 

a 

= 

a 

0 

+  kn+l 

Ja* 

2  v  n 

n  n+1  n+1 


n+1 


(1-  k2  )o2 
n+1  n 


These  are  the  celebrated  Levinson  recursions  for  sequentially  building  the 


whiteners  a,  or  equivalently  the  columns  of  the  A'  and  B'  in  the  Cholesky 
factorization  of  R~^. 

Vectorized  Fora  #1:  These  equations  nay  be  put  into  a  concise  matrix 
form  by  defining  the  following  matrices: 


J  A  J  =  Q  J  0  J  =  A 
Now  the  Levinson  recursions  are 
J-1  -  A  .»  ♦  0  J.” 

J.**1  -  H.1.  A 

A  cross  matrix  representation  for  these  recursions  is 

,n+l  m  r  A  mn 
*  Kn+1  A  i 

Jan+1  -  JKn+1  JJ  A  JJ  ia  -  Q  Jan 

where  the  cross  matrix  En+i  “  JKn+1J  is  defined  as  follows: 


Vhen  the  diagonals  cross,  the  center  term  is  (1+k  +j) 


This  representation  nay  be  used  over  and  over  again  to  obtain  the 


iterated  cross  aiatrix  recursion 

5°+1  “  *n+l  A  Kn  A  •  •  •  *2  A  Kj  A  1 


where  1  is  the  lxl  aiatrix  equal  to  unity 
inverted  to  give  a  handy  formula  for 
backwards: 


n  -1  n+1 

A  a  =  K  . ,  a 
n+1  - 


.  The  representation  may  also  be 
running  the  Levinson  recursions 


The  inverse  is 


-1  2  -1 
r  &  ( i  —  v  ) 

n+1  '  n+1' 


-k 


n+1 


-k 


n+1 


Of  course  k  *=  a°*J. 

n+1  n+1 


In  these  recursions,  a11  and  Jan  are  vectors  of 


increasing  dimension. 

Vectorized  Form  #2:  To  obtain  vectorized  equations  in  which  dimensions 
are  constant,  let  An  denote  a  t-vector  consisting  of  Jan,  followed  by  zeros, 
and  let  Bn  denote  a  t-vector  consisting  of  a11,  followed  by  zeros: 


The  vector  An  is  the  n**1  column  of  A',  but  Bn  is  not  the  n**1  column  of  B'. 
A”  *  JBn 


The  Levinson  recursions  produce  the  following  recsrsions  for  An  and  Bn: 


Bn+1  =  Bn  +  kn+1  T  An 
An+1  =  T  An  +  kn+1  B“ 

The  matrix  T  in  these  equations  is  a  delay  matrix: 


This  form  of  the  Levinson  recursions  will  simplify  our  discussion  of  fast 
algorithms  as  lattice  algorithms. 

Remarks :  1.  Given  the  reflection  coefficients  {k^  ^  1  ,  the  iterated 

cross-matrix  recursion  may  be  used  to  generate  all  of  the  order 
increasing  whiteners  an.  This  is  equivalent  to  factoring  R~*. 

2  t-1 

2.  Given  o.  ■  rn  and  the  reflection  coefficients  {k  J,  , 

u  u  ni 

2 

prediction  errors  a  may  be  generated  recursively. 

n 

3.  Given  the  highest  order  whitener  at+^>  with  kt+^  *  *t+l'  *^e 

cross-matrix  recursion  may  be  used  to  generate  all  of  the  order- 
decreasing  whiteners  a11,  and  all  of  the  order-decreasing  reflection 


coefficients  k  *  a  .  In  this  way  lattice  parameters  are  built  from 

A  A 


filter  coefficients.  The  stability  of  a  filter  built  from  a*  x  is 
tested  by  checking  for  |ka|  <  1  for  nsl»2, . . . , t-1 .  This  is  the  Shur- 


Cohn  test 


2 

4.  The  first  row  of  the  equation  BA'  =  HD  gives  the  following 
recursion  for  computing  correlations  from  whitening  filters: 
n 

2  .  T  n 
n  On  L  x  n- 1 
i=l 

These  remarks,  coupled  with  the  Levinson  recursions  themselves, 
complete  the  discussion  of  the  right  side  of  Figure  1. 

3.2  LeRoux-Gueguen  Algorithm 

Return  again  to  the  following  re-arrangement  of  the  UL  factorization  of  R: 
RA*  =  HD2 

Write  out  the  n**1  column  of  this  equation: 
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Ignore  the  zero  terms  in  the  right-hand  column  and  write 
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A  related  vector,  which  will  be  needed  in  our  development,  is  R^®: 
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The  first  tern  of  gn,  denoted  g&Q,  is  zero.  Note  the  definition  of  R^j  and 
summarize  these  results  as  follows: 
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Vectorized  Form  #1:  Our  objective  is  to  derive  recursions  for  h  ,  g  » 

2 

and  o  .  In  the  Appendix,  the  Levinson  recursions  for  Jan  are  used  to  derive 
n 

the  following  results: 
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Because  gn  #  Jh  ,  there  is  no  handy  cross-matrix  representation. 


Vectorized  Form  #2:  Non  define  Hn  to  be  a  t— vector  consisting  of 

2  n  2  n 

aD  5  >  preceded  by  zeros,  and  Gn  to  be  a  t-vector  consisting  of  g  . 

preceded  by  zeros: 


H“  =  o2 


n  h* 


Gn  -  c2  ' 

-  n  n 

i 


The  vector  Hn  is  the  ntfl  column  of  H,  but  Gn  is  not  the  ntE  column  of  G. 
Hn  ft  JGn. 

The  LeBoux-Gneguen  recursions  may  now  be  written 

gn+1  ,  gn  +  THn 


Hn+1  *  THn  +  kn+1  Gn 
where  T  is  still  the  delay  matrix. 


3.3  Summary 


In  order  to  emphasize  the  striking  similarity  in  the  Levinson-Durb in 
and  LeRoux-Gueguen  algorithms,  we  say  that  they  both  obey  the  recursion 

5n+i  -  U*  +  kn+1  TVn 
Vn+1  =  TVn  +  kn+1  Dn 

These  equations  are  reversed  (more  on  this  in  the  next  section)  as 
follows: 

gn  =  gn+1  _  k^+i  tv® 


TVn  +  k„+1  Gn 


The  iteration  for  increasing  order  filters  in  the  Levinson  recursions  is 
obtained  by  setting  initial  conditions  as  follows: 

V0  =  D0  =  (1.0,. ...0)' 

The  LeRouz-Gueguen  algorithm  is  obtained  by  setting 

20  =  *r0'rl» • • *»rt-l) '•  Yo  =  (0»ri» • • ' 

In  the  LeRouz-Gueguen  recursions,  the  variance  and  reflection 
coefficients  are  obtained  as  follows: 

=  (n+l)st  element  of  Va 
kn  =  -(n+2)nd  element  of  Un/cn^ 

These  recursions  are  discussed  as  vector  recursions  for  a  vector  processing 
machine  in  [10]. 


Figure  3  contains  four  different  lattice  cells:  (a)  and  (b)  are  un¬ 
normalized  end  normalized  forward  lattice  cells;  (c)  and  (d)  are  un¬ 


normalized  and  normalized  inverse  lattice  cells.  The  equations  implemented 
by  each  of  these  cells  are  illustrated  on  the  figure.  The  z_*  operator  is 
interpreted  as  the  delay  matrix  T  when  the  cells  map  vectors  Un  and  Vn  into 
vectors  Un+*  and  Vn+*.  When  the  vectors  Un  and  Vn  fill  up  sequentially  in 
time,  the  lattice  cells  see  scalar  sequences  and  the  z-^  operator  works  just 
like  a  scalar  delay.  When  these  cells  are  concatenated,  then  the  lattice 
structures  of  Figure  4  are  constructed. 

Denote  the  internal  variables  in  the  lattice  as  follows: 


‘t-lj 


V*  = 


'0 

n 


n 

lVi 


In  the  forward  (reverse)  lattice,  u°  is  the  internal  variable  at  the  output 

(input)  of  cell  n  at  time  i;  v°  is  the  input  to  cell  n+1  at  time  i  in  both 
latt ices . 

Think  of  the  forward  lattice  as  a  linear  transformation  that  maps  the 
input  sequence  U®  into  the  space-time  sequence 


(c)  inverse,  (d)  Variance-normalized  inverse 


while  generating  the  interval  variables  U  and  V  .  The  space-tine  sequence 
V  contains  upper  lattice  variable  at  cell  n  and  time  n  for  n=0, 1, . . . ,  t-1. 
Think  of  the  inverse  lattice  as  a  linear  transformation  that  maps  V  into  u® 
while  generating  the  internal  variables  Un  and  Vn. 


4.1  Forward  Lattice 


The  forward  lattice  may  be  used  to  (i)  convert  correlations  (r  }. 

n  u 

into  reflection  coefficients  {k^}^  variances  1»  end  "iaqmlse 


responses 


{h  q  *  in  the  factorization  I?  =HDH',  (ii)  convert  reflection 


coefficients  {k  *  into  variances  {o^}*  *  and  whiteners  {a*1}  in  the 
n  l  n  u  ~ 

-1  -2  t_l 

factorization  R  =  AD  A',  and  (iii)  whiten  correlated  data  {x  }n  to 

n  0 

obtain  uncorrelated  random  variables  in  the  analysis  representation  Du  =  Ax. 
The  Levinson  procedure  is  really  a  composite  procedure  consisting  of  (i)  and 
(ii)  above. 


( i)  Correlations  to  reflection  coefficients,  variances,  and  impulse 
re sponses .  Vhen  C,:  and  V®  are  set  with  the  initial  conditions  required  in 
the  LeRoux-Gueguen  algorithm  then  at  the  upper  and  lower  branches  of  the 
forward  lattice  one  observes 


Note  that  the  non-zero  entries  in  V  that  contain  information  about  h 
appear  as  trailing  entries  with  delay  n.  This  means  the  n**1  cell  need  not 

be  switched  into  the  lattice  until  time  n.  At  this  time  the  output  g”  *  in 


b“_1  -  G"_1 


n-1 


is  available  to  compute  kfi=  -g^  for  use  in  cell  n.  In 


summary,  the  correlation  sequence  (r^,  iy,  ...»  rt-l^  ®®y  t>®  fed 
sequentially  into  a  forward  lattice  in  which  cells  are  computed  and  switched 
in  sequentially  to  produce  reflection  coefficients,  variances,  and  *iaq>ulse 
2  n  n 

responses  o^h  .  The  elements  of  a  in  the  LeRoux-Gueguen  algorithm 

appear  sequentially  in  time  for  all  values  of  n  less  than  or  equal  to  the 
current  value  of  time.  Thus  the  lattice  procedure  produces  the  Cholesky 
factor  H,  row-by-row.  This  was  Friedlander ' s  insight.  The  procedure  is 
illustrated  in  Figure  5.  Note  that  initial  conditions  are  set  to  zero. 

To  complete  the  discussion,  we  note  that  one  can  stand  at  cell  n  and 

2  n 

observe  the  evolution  of  a  h  sequentially  in  time.  This  is  a  column-by¬ 


column  procedure. 

(ii)  Reflection  coefficients  to  variances  and  whiteners.  When  U®  =  V® 
=  (1 ,0, . . . ,0) ’ ,  then  at  the  upper  and  lower  branches  of  the  forward  lattice. 


one  observes 
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Figure  5 . 

Reading  Cholesky  Factors  out  of  the  lattice. 
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If  the  lattice  cells  are  wired  together,  as  in  Figure  4,  then  the  entries  in 
U®  =  (1,0,..., 0) '  may  be  sent  in  sequentially  and  the  entries  in  U®  and  V® 
may  be  read  out  sequentially.  This  follows  from  the  causal  dependence  of 
C®+*  and  Vn+1  on  C®  and  V®.  Thus,  given  the  reflection  coefficients  (kn) 
required  to  build  the  lattice,  an  impulse  may  be  sent  in  to  generate  order 
increasing  whiteners  Ja®  which  may  be  read  out  as  lead ing  entries  of 
V®  =  A®.  They  may  be  read  out  sequentially  to  produce  the  factorization  R~* 

A 

=  AD^A',  row-by-row.  But,  as  the  whiteners  a®  show  up  immediately  as 
leading  entries  of  A®  in  each  lattice  cell,  the  lattice  cannot  be 
sequentially  wired  as  in  Figure  5.  It  must  be  completely  wired  from  all  of 

the  reflection  coefficients  (k  }*  This  is  illustrated  in  Figure  6. 


An  alternative  is  to  generate  cells  sequentially,  and  then  after  each 
new  cell  is  computed  and  switched  into  the  lattice,  compute  the  impulse 
response  of  the  new  lattice.  Each  step  of  the  procedure  produces  A®  =  (Ja®, 
0)'.  This  composite  procedure  is,  in  fact,  a  Levinson-Durbin  procedure. 

(iii)  Correlated  data  to  uncorrelated  data.  If  the  forward  lattice  is 


« 


excited  with  correlated  inputs  (xq,xj 
in  U®  are 


zt-l^'  t^ea  variables  observed 


The  variable  u  is  the  n  entry  in  the  vector  Dn: 

n  - 

Du  =  Ax 

If  the  normalized  forward  lattice  is  excited  with  (xq, . . . »xt_j ) ,  then  the 
variable  u  n  is  the  entry  in  the  vector  n,  itself.  So,  by  reading  ont 


time  variables  u^,  u*,...,u^i  from  the  appropriate  cells,  at  the 

appropriate  times  (in  fact  the  time  equals  the  cell  nnmber),  we  read  ont  a 
white  sequence  of  nncorrelated  random  variables  with  normalized  variances. 
The  same  story  holds  for  the  nnnormalized  lattice,  except  now  the  variable 

n&  has  variance  o^.  This  Gram-Schmidt  orthogonal izat ion  is  illustrated 
n  n 

in  Figure  7. 

4.2  Inverse  Lattice 

The  inverse  lattice  may  be  used  to  (i)  convert  reflection  coefficients 

(k  }*  1  and  variances  {o^}*  1,  into  correlations  {r  *  and  'impulse 
n  u  n  u  n  u 

responses*  th11}*  1 ,  < i i )  convert  the  highest  order  whitener  a*  *  into 

reflection  coefficients  {k  * ,  and  lower  order  whiteners,  and  (iii)  color 

n  u 

uncorrelated  random  variables  Du  to  obtain  correlated  random  variables  in 
the  representation  x  *  HDu. 

(i)  Reflection  coefficients  and  variances  to  correlations  and  impulse 
responses.  From  Figure  5  we  see  that  initial  conditions  may  be  set  to  zero 
and  the  correlation  sequence  may  be  sent  into  the  lattice,  to  generate  the 
space-time  vector 


21 


This  means  initial  conditions  may  be  set  in  the  inverse  lattice  in  exactly 

2 

the  same  way  they  were  set  in  the  forward  lattice,  with  Oq  at  the  input  to 


the  first  cell  and  zeros  elsewhere,  to  produce  an  impulse  response  for 


the  inverse  lattice  which  generates  the  internal  variables  Ha  -  (0,  hn) 

and  the  ontpat  seqnence  {r^}*  *.  This  generalizes  the  Robinson-Treitel 

result,  to  include  the  generation  of  Cholesky  factors  as  internal  lattice 
variables.  It  is  illustrated  in  Figure  8. 

(ii)  Highest  order  whitener  to  reflection  coefficients  and  lower  order 
whiteners.  From  Figure  6  we  see  that  initial  conditions  may  be  set  to  zero, 
and  a  1  placed  at  the  input  to  the  first  cell,  to  produce  the  order- 
increasing  whiteners  tn  (or  Jan)  as  internal  variables  in  the  lattice.  The 
ouput  of  the  (t-l)st  cell  contains  Bt_*  =  at_*  on  the  upper  branch.  This 
means  the  initial  conditions  may  be  set  to  zero,  and  the  sequence 

(1,  a*  *,  ...,  a**)  sent  into  the  inverse  lattice  to  produce  all  of  the 

lower  order  whiteners  as  internal  variables.  The  lattice  must  be  fnlly 
connected  to  achieve  this.  See  Figure  9. 

(iii)  Dncorrelated  data  to  correlated  data.  Recall  Figure  7: 
with  initial  conditions  set  to  zero  and  the  input  equal  to 

(xft,  x,,...,x^  , ),  the  uncorrelated  random  variables  u11  are  observed  at  cell 

n  at  time  n.  Thus,  an  inverse  lattice  may  be  set  with  zero  initial 

conditions  and  excited  with  uncorrelated  random  variables  u~.  uj, . . .  .u11, . . . , 

u  l  n 

u*_j  ,  each  entered  at  the  appropriate  time  and  place  (cell  n,  time  n),  to 

produce  the  correlated  output  (xq,  xj,  . . .  .x^j) .  This  is  illustrated  in 
Figure  10. 
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Figure  8 . 

Generating  Cholesky  Factors  and  Correlations  from 
Reflection  Coefficients. 


5.  Autoregressive  Moving  Average  Sequences 
So  far  we  have  dealt  exclusively  with  a  wide-sense  stationary  time 
series  (xn)  whose  correlation  sequence  {rnJ  is  even  and  non-negative 

definite.  A  finite  record  of  the  series  {x}  „  *  has  a  correlation  matrix 

n  0 

R  which  is  symmetric,  Toeplitz,  and  non-negative  definite.  We  now 
specialize  our  results  to  the  case  where  the  correlation  sequence  obeys  the 
coupled  autoregressive  moving  average  (ARMA)  recursions 


)  a  r  =  >  b  h  ;  a.  =  1 

4  n  t-n  4  n  t+n  0 


3  a  h  =  b 

4  n  t-n  t 


s  b0  =  b0 


Such  a  correlation  sequence  is  said  to  be  ARMA  (p,p). 

Whenever  these  recursions  hold  for  the  correlation  sequence,  the 
underlying  time  series  itself  is  said  to  be  ARMA  (p,p).  It  obeys  the 
autoregressive  moving  average  recursion 


^  a  x  =  ^  b  u 

4  n  t-n  4  n  t-i 


{uj.}:  sequence  of  uncorrelation  random  variables  with  unit  variance. 

The  operator  or  transfer  function  representation  of  {xnJ  is 
A(z) {x.}  *  B(z)  {u*} 


p  p 

A(z)  =  ^  a^z  B(z)  =  ^  b^z  n 

n=0  n=0 

This  may  also  be  written  as  an  infinite  moving  average, 
{xtJ  =  H(z)  {nt)  , 

where  the  transfer  function  H(z)  obeys  the  recursion 
A(z )H(z)  =  B(z) 

P 

5  a  h  =  b 
Z  n  t-n  t 

n=0 

The  spectrum  E(z)  =  H(z)H(z-*)  obeys  the  recursion 


A(z)R(z)  =  B(z)H(z  1) 

P  P 

\  a  r  =  5  b  h 

Z  a  t-n  Z  n  t+n 

n=0  n=0 

5.1  Stationary  Wold  Representation 

The  sequence  {xnJ  may  be  decomposed  as  follows: 

xt  =  xt/t-l  +  Vt 

xt/t-l  =  (H(z)  '  h0)nt 

=  (B(z)/A(z)  -  hg}ut 

This  may  be  rewritten  as 

A(z )  =  (B(z)  -  A(z)h0)  (utJ 

xt  =  xt/t-l  +  h0  ut 
Define  the  polynomials 

z~^Q(z)  =  B(z)  -  ho  A(z) 
z"1  P(z)  =  1  -  A(z) 

The  previous  equations  can  then  be  re-cast  as  follows: 
xt+l/t  =  P(z)  xt/t-l  +  Q(z)ot 


c 


xt  "  xt/t-l  +  h0  Bt 


A  block  diagram  is  illustrated  in  Figure  11.  The  variable  is  the 

minimum  mean-squared  error  prediction  of  xt  based  on  the  infinite  past 

{ • • . U^ ,  • • • i  U  ^  • 

For  this  structure  to  be  used  as  a  synthesizer  of  a  stationary  time 
series,  the  inputs  u^  must  be  initiated  infinitely  far  in  the  past. 
Alternatively,  the  initial  conditions  in  the  structure  must  be  set 
appropriately.  This  brings  us  to  the  Markovian  representation. 


5.2  Markovian  Representation 

The  Markovian  state  space  representation  corresponding  to  the 
stationary  Void  representation  is 

it+1  =  F  xt  +  ^  Bt 
*t  =  §'  xt  +  h0  ut 

ut:  sequence  of  uncorrelated,  unit  variance  random  variables 


In  this  state  space  model,  the  vectors  and  matrices  are  defined  as  follows: 
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t+p+i/t-i. 
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h  with  a  time-varying  Salman  gain  k^,  replace  the  random  initial  conditions 
with  zero  initial  conditions,  and  replace  the  nnit  variance  inpnt  sequence 
{Ujj}  with  a  sequence  {onun}  whose  variance  is  time  varying.  The 
representation  is 

?t+l  =  F»t  +  *t  °t  nt 

xt  =  it  +  h0  fft  nt 

ut:  sequence  of  uncorrelated  random  variables  with  unit  variance 

This  representation  is  illustrated  in  Figure  13. 

For  the  innovations  representation  to  produce  a  stationary  sequence,  we 
require  the  correlation  of  {z^}  to  be  {rQ}.  That  is,  we  require 


The  expression  for  r*  is 

r*  =  6'  Fn  Qt  §  +  hQ  6'  Fn_1  kt 

where  Qt  is  the  zero-lag  state  covariance  at  time  t,  for  this  to  equal  rt  we 
require 

kt  ot  =  -F(Qt-QQ)  +  h 


A  -  «’<W*  +  A  -  *0  -  >■  i 

This  is  the  usual  Ricatti  solution  given  for  the  Kalman  gain  k( 

2 

and  innovations  variance  o .  But  there  is  another  way  to  go. 

In  the  innovations  representation,  the  output  may  be  written  as  a  time- 
varying  convolution  of  a  time-varying  impulse  response  with  the  input  (un) : 


The  impulse  response  is 

0  ,  t  <  n 

hl-n  =  h0  “  1  '  1  =  « 

6'  Ft_n_1  k  ,  t  >  n 
— 

This  means  the  vector  of  outputs  x  =  (xq,  Xj,  xt_i> *  may  be  written 


or  as  follows: 


So,  in  fact,  the  time  varying  impulse  response  in  the  innovations 
representation  just  fills  in  the  Cholesky  factorisation 
R  -  HD2H' 

corresponding  to  the  synthesis 
x  -  HDu 

Why  is  this  important?  Remember  we  have  a  fast  algorithm  for  obtaining  the 


columns  (or  rows)  of  H.  And,  from  the  expression  that  relates  the  time 
varying  impulse  response  to  the  Kalman  gains  we  have  the  resnlt 


This  means  we  can  rnn  a  fast  Cholesky  algorithm  (either  by  columns  or  rows) 

and  read  ont  the  Kalman  gains  k  as  the  entries  h°  through  h&.  This  is  a 

-n  1  p 

fast  Kalman  gain  algorithm.  An  attractive  iaplementation  uses  the  lattice 

excited  with  a  correlation  sequence  {r^}*  1.  The  gains  are  read  out  as 

internal  lattice  variables.  This  is  illustrated  in  Figure  14.  Of  course 
the  gains  kt  inherit  the  same  recursions  as  hn  and  gn  in  Chapter  4. 

As  the  final  topper  to  this  story,  note  that  the  companion  matrix  F 
satisfies  its  own  characteristic  equation: 
n+p 

}  a  Fn+P_t-  o 
L  t-n 


Thus  the  time  varying  impulse  response  h”_n  obeys  the  following  recursion: 


n+p  n+p 

»-  }  .  «  3  . 

L  t-n  -n  L  t- 

t**n  t=n 


n  n+p-t+1 


n+p+1 

■  \  •, 


L  p+l-t+n  h° 

.  t-n 

t=n+l 


This  means  each  coluam  hn  may  be  generated  by  cooputing  the  variance  o^,  the 


Kalman  gain  k^,  and  then  using  the  gain  to  initialize  the  recursion  above  to 

fill  out  the  columns.  The  same  recursions  can  be  derived  for  Vi  «’)• 

leads  to  the  Morf,  Sidhu,  Kailath  recursions  for  k  [9]. 

-  n 

6.  Conclusions 

Our  conclusions  are  much  like  our  introductory  comments.  Forward  and 
inverse  lattices  may  be  used  with  a  variety  of  excitations  and  initial 
conditions  to  produce  internal  lattice  variables  which  are  Cholesky  factors 
of  correlation  matrices  and  their  inverses.  When  the  correlation  matrix  is 
ARMA,  then  these  internal  variables  may  be  used  to  generate  Kalman  gains  and 
to  identify  autoregressive  parameters.  The  highest  order  Kalman  gain  may  be 
identified  with  a  stationary  impulse  response  to  obtain  an  estimate  of  the 
moving  average  coefficients. 

In  references  [10]  the  vector  recursions  of  Section  3.0  are  presented 
as  algorithms  that  may  be  implemented  on  a  vector  processing  machine.  In 
[11]  a  fixed  length,  time  varying  lattice  is  derived  for  implementing  the 
Morf,  Sidhu,  Kailath  recursions. 


8.  Appendix 

Begin  with  the  UL  factorization  of  R: 

RA’  =  HD2 

Write  out  the  n**1  column  of  this  equation: 


r, .  .  .  r 

0  1  n 

rl  r2*  •  *  rn+l 

Ja“ 

2 

0 

=o 

•  • 

•  • 

n 

•  • 

rt-l  ‘  •  •  rt-l-n. 

h” 

Ignore  the  zero  terms  in  the  right-hand  column  and  write 


A  related  vector,  is  R_an: 


The  first  term  of  gn,  denoted  SDq.  is  zero.  Note  the  definition  of  and 

summarize  these  results  as  follows: 


n  n  2 

Our  objective  is  to  derive  recursions  for  h  ,  g  ,  and  o  .  To  this  end, 

—  n 


2  n+1 

add  one  more  row  to  the  equation  governing  £  : 


2  n+1 

Vi  *t 


n+1  n+1  , 

rt  a  +  r  a  +  . 
t  n+1  t-1  n 


,  n+1 

rt-n-l  *0 


Use  the  Levinson-Dnrbin  recursion  for  Ja 


From  this  recursion  we  note  the  following  recursion  for  a" 


2  2  ,,  ai.  2.,  ,  n. 

o  =  o  (1-k  g«)  =  o  (1-k  , ,  g.) 

n+1  n  n+1  *0  n  n+1  “1 


To  obtain  a  recursion  for  gn,  consider 


Add  one  more  row  and  use  the  Levinson-Durbin  recursions  to  obtain 


From  this  recursion  we  note 


n+1 


An  n 
~*0  =  _gl 


Substitute  into  the  recursion  for  a  to  obtain 

n 


n+1 


o2  (1-k2  ) 

n  n+1 


Ve  may  summarize  the  recursions  as  follows: 


2 

0 

.n+1 

2 

0  ‘ 

,  n 

^  .  2 

n 

°n+l 

h 

=  a 

n 

h 

+  h  o 
n+1  n 

£ 

>rl. 

• 

n.  2 
c.  /o 

t  nj 

2 

0  • 

n+1 

2 

n 

.  2 

0  ‘ 

.  n 

°n+l 

£ 

=  a 

n 

£ 

+  k  ,o 
n+1  n 

h 

n+1 

K 

n 

.  "t 

If  we  keep  only  the  relevant  terms 


and  Q  we 

may  write 

2 

.  n+1 

2 

°n+l 

b  =  o 

n 

2 

n+1 

2 

ffn+l 

£ 

=  o 

n 

A'  hD  +  k  -  a2 
n+1  n 


k  , -  a 
n+1  n 


>  and  use  our  previous  definitions  of 


Q'  £n 
A'  hn 


These  recursions  are  initialized  as  follows: 


r. 

r 

,  • 


t--, 

L "  . 

T 


2  / 

°o  5  "  <ro* 


2  0 


rt_i)'  •  aQ  t  =  (0, r1# . . . , rtl) ' 
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